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ABSTRACT 


This  paper  presents  a  method  of  optimal  filter 
design  for  sampled  data  systems,  based  on  the  theory 
originally  developed  by  R.  E»  Kalman.  The  first  half 
of  the  paper  deals  with  the  theoretical  development 
of  mathematical  models  for  linear,  discrete  dynamic 
processes  and  the  optimal  filter  equations  for  such 
processes.  The  latter  half  discusses  digital  pro¬ 
gramming  techniques  for  optimal  filter  design  followed 
by  two  illustrative  examples. 


PREFACE 


During  the  past  several  decades  a  fair  amount  of 
theoretical  effort  has  been  devoted  to  the  study  of 
problems  ■which  are  of  a  statistical  nature.  Not  the 
least  important  is  the  class  of  problems  in  communica¬ 
tion  and  control  which  involves  the  separation  of  random 
signals  from  random  noise,  or  the  estimation  of  the 
states  of  a  dynamic  process  based  on  noisy  observations 
of  a  few  of  these  states. 

In  several  papers  written  since  I960,  R.  E. 
Kalman  developed  a  theoretical  approach  for  optimization 
of  filters  for  the  above  mentioned  class  of  problems. 
The  theory  is  not  all-embracing  In  that  certain  con¬ 
ditions  must  be  satisfied  before  his  technique  can  be 
employed. 

The  intention  of  this  paper  is  to  present  a  method 
for  optimal  filter  design  for  sampled  data  systems, 
based  on  Kalman’s  approach.  The  first  half  of  the 
paper  deals  with  the  theoretical  development  of  mathe¬ 
matical  model  parameters  for  linear  dynamic  processes 


and  the  optimal  filter  equations  for  such  processes 


The  latter  half  discusses  digital  programming  techniques 
for  optimal  filter  design  followed  by  two  illustrative 
examples . 
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CHAPTER  I 


MODELS  FOR  RANDOM  PROCESSES 

Before  one  can  hope  to  achieve  any  amount  of 
effective  filtering,  it  is  necessary  that  a  fair  amount  of 
knowledge,  about  the  physical  phenomena  to  be  observed, 
is  known.  For  instance,  if  a  sine  wave  buried  in  noise, 
is  to  be  recovered,  an  apriori  knowledge  of  the  signal, 
i.e.  frequency  of  the  sine  wave,  is  necessary.  In  addi¬ 
tion  if  the  statistics  of  the  noise  are  known  then 
optimum  filtering  can  be  achieved.  It  therefore  becomes 
necessary  to  make  a  study  of  the  message  (signal) 
process  before  the  construction  of  a  filter  is  attempted. 
To  maintain  generality  we  will  henceforth  only  consider 
random  signals  with  the  added  stipulation  that  these 
signals  are  produced  by  linear  dynamic  systems  excited 
by  white  noise. 

1-1  LINEAR  DYNAMIC' SYSTEMS  (CONTINUOUS 
TIME)  . 

Since  we  are  concerned  only  with  linear  dynamic 
systems  a  brief  review  of  linear  differential  equations  is 


in  order. 


A  first  order  differential  equation 


+  C(  X  =  U 


(1.1) 


has  a  solution  (see  Appendix  I) 

-e^Zo  -t  } e^(t^ u<r)Jr  {1,2) 

0 

where  £  is  the  homogenous  solution  and 

o 

is  the  particular  solution. 

Consider  now  a  set  of  first  order  differential 
equations,  which  define  a  linear  dynamical  systems 

^  «)  <1,3) 


or  In  vector  notation 


JL  s  f(i)X  +  D(t)  u(t) 


(1.4) 


where  x  and  u  are  1  x  n  column  vectors  and  F  and  D 
are  n  x  n  matrices. 

The  solution  (see  Appendix  I)  to  this  set  of 
equations  is: 

xw  =  *Xo  h»  cir 

0 


or  it  may  be  written 

t 

+  J  §  Lt,r)0(r)  u(r)  dr  (1-6> 

By  definition  we  call  the  vector  x  the  state  of 
the  system  and  u  the  input  or  control  function. 

Since  all  states  (x-L)  may  not  be  observable  we 
define  the  output  of  the  system  to  be 

at*)  -  hW  *  <y  (K7) 

where  ^(t)  is  a  p  vector  and  H(t)  is  a  pxn  matrix. 
If  all  states  were  observable  then  H  would  be  equal  to 
the  identity  matrix  I . 


We  can  now  represent  the  system  in  matrix  block 
diagram  form  as  shown  in  Pig.  1-1.  ' 


Fig.  1-1  Matrix  block  diagram  of  a  linear  dynamic  system. 


The  integrator  in  Fig.  1-1  actually  represents  n 
integrators,  one  for  each  state  of  the  system,  while 
F(t)  shows  how  the  outputs  of  the  integrators  are  fed 
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back  to  the  inputs  of  the  various  integrators.  Perhaps 
a  look  at  a  simple  2-state  system  at  this  time  might 
clarify  Fig.  1-1. 

Given  the  linear  dynamic  system  of  Fig.  1-2,  deter¬ 
mine  F{  t ) ,  D(t),  and  H(t). 


Fig.  1-2  A  2-state  system 


We  can  immediately  write  down  the  equations  for 
the  system  as: 


X,(t) 

=  xi  (V 

(1.8) 

-  -  -  o(  XL{tJ  f  u.(tj 

(1.9) 

and  our 

observable  state(s) 

y  w 

(1.10) 

It  is  immediately  obvious  that 


thus  giving  the  vector  differential  equation 


k 


and 


a  =  L‘  °J  \x' 

Lx* 


(1.12) 


1-2  LINEAR  DYNAMIC  SYSTL-  '  •C'R 
TIME)  . 

If  we  specify  the  equations  of  a  linear  dynamic 
system  in  the  form  of  difference  equations  then  they 
are  easily  mechanized  on  present  day  digital  computers. 
With  this  in  mind  the  scope  of  this  paper  will  be 
directed  towards  discrete-time  situations. 

In  (1.6)  we  see  that  the  continuous-time  solution 
to  a  linear  dynamic  system  is: 

xCO  =  x(t,)  +  J  (j)(t,rj  D(rJ  g(r)dT  0.6) 

If  u(r)  is  held  constant  over  the  interval  of 
integration  then  we  obtain: 

x(t)  z  $(t,t0)£fa)  +  A  (t  ,-t0)  uW  (1.13) 

where 

,*  . 

=  J  $  (t,r)  D(r)  dr  (1.14) 

•to 

or  more  conveniently 

x(tfi)  =  (p  (t  th  t)  iW  +  A  (1.15) 
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In  (1.15)  we  assume  a  sampling  period  of  one  time 
unit.  A  block  diagram  of  the  linear  discrete-dynamic 
system  is  shown  in  Fig.  1-3 • 


Fig.  1-3  Matrix  block  diagram  of  a  linear  discrete- 
dynamic  system. 

1-3  DETERMINATION  OF  MODEL  PARAMETERS. 

The  matrix  <!>(**<*  t)  oc curing  in  (1.6)  and  (1.13)  is 
called  the  TRANSITION  matrix  and  has  the  following 


properties: 

4?  (to,  to)  U  I  (Identity  Matrix) 

(1.16) 

^  (tz  >  t})  d)  (t  1  y  to)  -  Q  (^2  )  to) 

(1.17) 

=  F(t) 

(1.18) 

at 

(1.16)  and  (1.17)  are  fairly  obvious  and  (1.18)  is 
obtained  by  setting  u(t)  to  zero  and  differentiating 
(1.6).  These  properties  can  be  useful  in  checking  the 
accuracy  of  analytic  expressions  for  the  $  matrix. 

If  the  F  .matrix  is  constant  then  the  transition 
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matrix  elements  depend  only  on  the  time  difference 
t-t0  and  can  be  calculated  from  the  following  expres¬ 
sion: 

$  (t>t.)  =  e*  -  £  [F(t-t.)j L/i\  (1.19) 

L-o 

A  second  (and  easier)  method  for  obtaining  $(t,t0) 
is  by  the  use  of  signal  flow  techniques  and  the  appli¬ 
cation  of  unit  impulses  to  the  input  of  selected 
integrators.  This  method  will  be  demonstrated  else¬ 
where  in  this  paper. 

The  A  matrix  may  be  obtained  by  performing  the 
integration  in  (1.11+)  or  by  the  second  method  mentioned 
above  for  the  (£>  matrix. 

1-1+  THE  GAUSS-MARKOV  PROPERTY. 

A  large  number  of  physical  phenomena  possess  the 
Markov  property.  Basically  it  means  that  the  best 
estimate  of  a  future  state  of  a  process  can  be  made 
without  the  knowledge  of  all  the  past  history  of  the 
process.  In  a  strict  sense  it  implies  that  the  best 
estimate  of  a  future  state  can  be  derived  from  the 
last  observation  of  the  states.  A  very  trivial  example 
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•would  be  the  motion  of  a  particle  with  a  constant  velocity 
vector.  Given  the  best  estimate  of  the  present  position 
and  velocity  of  the  particle  one  can  formulate  a  best 
estimate  of  position  and  velocity  for  any  time  in  the 
future.  In  fact  the  output  from  any  linear  dynamic 
system  is  Markovian.  If  u(t)  Is  set  equal  to  zero  in 
(1.15)  then  this  property  may  be  expressed  mathematically 
as: 

X  (t+0  =<$  (1.20) 

If  u(t)  is  a  gaussian  random  vector  then  the  sequence 
of  random  vectors  ...x(t-l),  x(t/l),  ...  generated  by 
(1.15)  is  known  as  a  gausc-Markov  sequence.  The  stipu¬ 
lation  that  u(t)  is  gaussian  implies  that  the  sequence  ...r 
u(t-l),  u(t),  u(t/l),  ...  are  normally  distributed  random 
vectors  such  that  the  cross-variance  matrix: 

COV  [u(t,)  J  u  (te)J  Z  o  iort,ttz  (1.21) 

I.e.  u(ti)  and  u(t2)  are  independent.  In  addition  the 
random  vectors  are  completely  defined  by  specifying  their 
first  and  second  order  moments,  i.e.  E(u(t))  and 
B(u(t).u(t)  ).  For  the  purposes  of  this  paper  u(t) 
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will  be  assumed  to  have  zero  mean; 

i.  e.  E  [  U  ( t  )j  =•  0  for  all  t  (1.32) 

E(u(t).u(t)  )  is  called  the  auto-covariance  matrix 
of  the  vector  u(t)  and  will  be  denoted  by  U(t),  i.e: 

si  =  u (?)  d-23) 

or 

cov  [  u  (t)  j  t  U(t) 

Considering  now  the  state  of  a  process,  we  assume 
that  the  initial  state,  x(t0),  is  a  gaussian  random 
variable  of  zero  mean  and  arbitrary  variance.  By 
repeated  application  of  (1.15) »  we  see  that  future  states, 
. ..,  x(t-l),  x(t),  x(t/l),  ...  will  also  be  gaussian 
random  variables,  since  they  are  obtained  by  linear  com¬ 
binations  of  gaussian  random  variables. 

In  probability  terminology  we  may  now  define  the 
Gauss-Markov  property.  Since  u(ti)  and  u(tg)  are 
independent  for  tj  yfc  tg»  then  the  conditional  probability 
distribution  of  x(t)  is  dependent  only  on  the  previous 
state,  i.e: 

p(£(t)  £  ll|2C(t-l)JX  (t-a),X  (t-3); - )  (1.24) 

n  P  (x(t)  <  |  X(t-I)) 
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Where  £j_  is  an  arbitrary  vector. 

1-5  THE  COMPLETE  MODEL 

In  Fig.  1-3  the  observable  output  vector  ^(t)  cannot 
be  measured  with  infinite  accuracy  and  therefore  to  com¬ 
plete  the  model  for  random  processes  (with  previously 
mentioned  restrictions)  ,  a  source  of  measurement  noise 
must  be  added.  This  is  illustrated  in  Fig.  1-2+  where  v(t) 


Fig.  1-2+  Model  for  random  processes  generated  by  discrete¬ 
time  linear  dynamic  systems. 

as  ^(t).  v(t)  *s  w^te  noise  (gaussian),  which  we  assume 
to  have  zero  mean  with  arbitrary  variance: 

£  [  V*(t)]  -  0  (1.25) 

£  GO  =  C0\j[_vit)\  =  R(t)  (1.26) 

In  addition  we  specify  that  v(ty)  and  v(t2)  are  indepen¬ 
dent  for  1 1  y  t2 1  i.e: 

COV[v(t0;l/(ti)]  =0  for  t.*tt  (1.27) 
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CHAPTER  II 


THE  KALMAN  FILTER 

2-1  DEFINITION  OF  THE  FILTERING  PROBLEM 

In  Chapter  I,  a  model  of  a  linear  dynamic  system, 
excited  by  white  noise,  was  developed.  The  purpose  of 
the  Kalman  filter  is  to  give  a  best  estimate  of  all  states 
of  the  system,  based  on  noisy  observations  of  the 
observable  states.  Since  the  system  is  linear  we  may 
write 

t)  =  X  <0  +  GCt)  U  (t)  -£(*)]  (2.1) 

where  x*(t)  is  the  best  estimate  of  x(t),  based  on  the 

\ 

current  observation  z^t), 

x(t)  is  the  best  estimate  of  x(t),  based  on  the 
previous  observation  z^(  t-1) , 

£(t)  is  the  best  estimate  of  z(t),  based  on  the 
previous  observation  z_(t— 1),  and, 

G(t)  is  an  nxp  gain  matrix,  the  magnitude  of  its 

< 

elements  being  indicative  of  the  amount  of  information 
carried  in  z(  t)  . 
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Tr  solve  the  filtering  problem,  the  filter  must 


therefore  determine  values  for  the  three  unknowns  on 
the  right  hand  side  of  (2-1),  namely  x(t),  _z(t),  and  G(t). 
2-2  SOLUTION  OF  THE  FILTERING  PROBLEM. 

Since  we  assume  a  complete  knowledge  of  the 
dynamics  of  the  system,  computation  of  x(t)  and  is 

quite  simple. 

&(t)  -  e  [*(*)! 

=  (t.t-l)  +  (2-2) 

however  since  0  for  all  t  then 

X(c)  =  (2*3) 

In  the  model  we  saw  that 

g(t)  =  %  it)  +  of  CO  (1.28) 

r  HC+)  x(t)  +  \[(t) 

now  2(t)  Z  ££g(^J  |g(t~/)J 

=  Hit)  I  2(t-0]  +££  iC(D] 

=  H  (t)  X.  (t)  (2.4) 

since  £  [if  CP)]  =  & 
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Having  now  developed  expressions  for  x(t)and  z(t) 
a  matrix  block  diagram  of  the  filter  (see  Pig.  2-1)  can 
be  produced.  The  only  unknown  yet  to  be  calculated  is 
the  optimal  gain  matrix  G(t).  Before  approaching  this 
'calculation  a  criterion  for  optimal  must  be  specified. 


Fig.  2-1  THE  OPTIMAL  FILTER. 

The  criterion  used  is  that  we  wish  to  find  G(t) 


such  that  the  loss  function 

L  z  E  [CXCt)  -X*ft))T(X (t) - ^  *))J  is  minimized.  (2.5) 
That  is  to  s'ay  that  the  sum  of  the  variances  of  the 
errors  associated  with  the  estimate  of  the  individual 
states  is  minimized.  Because  the  errors  are  gaussian 
it  can  be  shown  (ref.  1)  that  this  criterion  will  in  faot 
produoe  an  optimal  gain  matrix. 

A  number  of  different  derivations  for  G(t)  are 
available  in  the  literature.  For  the  most  part  these 
derivations  are  mathematically  rigorous  and  somewhat 
complex.  Perhaps  the  easiest  one  to  follow  is  a  semi- 
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We  wish  to 


heuristic  approach  used  by  Schmidt, 
minimise  the  scalar  loss  function  defined  by  (2.5)  •  Note 


E[£c  -  >c*)r  ( Z'X?)] --  Trace  eux-^xx-Z*)]  <2-6) 


where  the  Trace  operator  denotes  the  sum  of  the  main 
diagonal  elements,  and  (t)  is  left  out  to  avoid  unneces¬ 
sary  clutter. 

Expanding  the  covariance  matrix  in  (2.6),  using 
(2.1)  ,  We  obtain 

Efe-x*>tl£-Z*) f)=  BlU&-x)-S(g-i)]l(M-S)-&(l-i)]  ] 

-  E[<x-i)(x-£)Tj  -  E[n(*-i)te-£)rJ 

-£[pc  4)(z-i)V]  +  e[6  fe-iXs-H  )VJ 

but  £  :  tlT  )  >  £  r  H  £ 

Substituting  for  z  and  _£,  and  noting  that  £[(X~X)ifT]*0 
we  obtain 

«  Lr[<X-2)Cx4)T]-E[GR  Qe-tj&S/] 
*-eL(X-&)Uim£)rHG] 
ciG  ( H  (x x  )Iht*  ¥‘!£r)Gr] 
-P-CHP  -  PHTGT  t  G  (H  PHT+  ft)Gr  (2.7) 

where  Ps  £  [(X'm£)(X-£)r J  and  =  £["if 'if  r_J 


We  now  wish  to  find  an  expression  for  G  such  that 
the  trace  of  (2.7)  is  a  minimum.  Since  the  terms  In 
(2.7)  are  matrices  this  could  become  a  very  arduous  task. 
Let  us  consider  for  a  moment  that  (2.7)  Is  a  scalar 
expression  (i.e:  the  matrices  are  lxl  in  dimension),  thus 
reducing  the  right  hand  side  to: 

P  -Z6PHr+  6*  (H  PHt  +  R) 

We  now  differentiate  with  respect  to  G  and  set  the 
resultant  to  zero  obtaining 

-  2  PHT  +  2  (HPHt+  R)G  =  0 
or  S  s  P  Bt  (H  PHt  +  R )*'  (2.8) 

It  can  now  be  shown  that  (2.8)  will  in  general 
provide  the  optimum  gain  matrix  by  letting 
C  =  G  -  PWT(HPHTtR)~' 

or  6  -  C  +  P  Hr  C  H  P  h7  +  ft)'1  (2.9) 


Simple  substitution  of  (2.9)  for  G  In  (2.7)  will 
reveal  that  the  trace  of  (2.7)  will  be  a  minimum  for 
Cb*Q.  Thun  (2.8)  provides  us  with  the  optimum  gain 
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equation* 


Combining  equations  (2.7)  and  (2.8)  results  In  an 
expression  for  the  covariance  matrix  of  the  error  in  the 
filter’s  estimate  of  the  states  of  the  system; 

£[(X-**)(X-X*)TJs  P-PHt(H  PHT+R)~iH P  -  PHT  [PHT  (  H  PHrtR)'‘]T 

'  f£pHT(HPHT  +  «r/J[HPHT+Rj[PHVHPH+R)"'JT 

=  Pit) - p(t)nh Gwo m  (2.10) 

^  f  -  Q?  c  (jjr  &  f* 

In  order  to  complete  the  filtering  problem  a  recur¬ 


sive  relation  for  the  conditional  covariance  matrix, 
P(t,t-l)jr  must  be  derived. 


Recall  that 

but  z(4)  t 

and  a  d)  (t)t-l)  X? (tmi 

/.  P(t,t-|)  =  E{[#(^ t-0(a Lt-i)- -**(*-/,*-/))  +  Mt>t -On M]*  y 


(2.11) 

+  ae[^tJ 


rlAT 


?  b 


/j-oJr  A>f  (W.x^c  ~"A 


o-a  cl  ^-uaaa-Ca^a. 


4-  /i-  * 


Bquation  (2.13)  is  called  the  variance  equation  which  . 
denotes  the  covariance  of  the  error  between  the  actual 
states  x(t),  and  the  predicted  states  x(t,t-l). 

Since  (2.13).  is  recursive,  an  initial  covariance  matrix, 
P(tQ),  must  be  specified,  and,  since  we  assume  that  x(tG) 
is  gaussian,  with  zero  mean,  our  best  estimate  of  x(tQ)  is 
zero.  Hence 

P(t.)  e  f  [x  (2.14) 


In  determining  F(t0),  one  will  often  find  that  the 
elements  off.  the  main  diagonal  will  be  zero,  that  is  to 
say  that  the  individual  initial  states  are  independent  of 
one  another.  To  illustrate  this,  perhaps  a  simple  example 
will  be  helpful.  Let  us  suppose  that  we  are  going  to  make 
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observations  of  the  position  of  a  particle  in  motion,  and 
assume  that  the  particle's  velocity  is  constant  but  unknown* 
The  two  states  of  the  system  are  position  and  velocity. 

Any  knowledge  one  might  have  about  one  of  the  initial 
states  will  in  no  -way  assist  in  determining  the  other  initial 
state.  Hence  the  two  states  are  initially  independent  inso¬ 
far  as  the  observer  is  concerned.  However,  as  more 
observations  are  made,  the  two  states  build  up  a  depend¬ 
ency  and  the  off-diagonal  elements  in  P(t)  ,  generated  by 
(2.13),  become  non-zero. 

2-3.  SUMMARY  OF  FILTER  EQUATIONS. 

.  For  convenience  the  equations  for  the  optimal  filter 
are  grouped  belpws  J 

Vh^.U)  =i>C£,  +  9<t>  1  ^ 

\ Gft)  =  P(J)  Hr(t)  [h  (t)P(t)(jr(t)  +  RCt )]~*  II  | 

ii.ftjt-i)  =  §  (t)t-i)  '  hi  } 

t-0  =  H(t)  t-/)  v  iv  1 

)  =  X  (tjt-l)  +  G(t)  [z6tr)-j[(tjt-/ji]  v 


'  0-1- 
c  ;  J 


Ct-v . vA  •  A.  v  A  • 

'  ■  0 


uifl  CmM 


S-  P  f+  r  j 
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CHAPTER  III 


FILTER  DESIGN 

3-1  SOME  PRELIMINARY  CONSIDERATIONS . 

The  equations  for  the  optimal  filter  were  derived  in 
Chapter  II,  and  are  summarized  In  Section  2-3  •  Design 
of  the  optimal  filter  consists  mainly  of  writing  a  suitable 
computer  program  to  carry  out  the  calculations  indicated 
by  the  filter  equations  I  through  V.  The  Input  to  the 
filter  is  the  noisy  observation  vector,  z.(t),  and  the  out¬ 
put  is  the  best  estimate  of  the  system  state  vector, 
x*(t). 

One  of  the  chief  problems  In  carring  out  the  filter 
computations  is  the  determination  of  the  inverse  of  the 
matrix  found  in  equation  II.  If  this  matrix  Is  singular, 
then  the  inverse  does  not  exist.  One  must  then  resort 
to  the  calculation  of  a  pseudo-inverse.  The  manner  in 
which  the  pseudo-inverse  is  determined  is  shown  in  [2] . 
and  will  not  be  discussed  further  here.  In  either  case 
the  time  required  for  Inverse  computation  is  relatively 
high,  the  time  being  roughly  proportional  to  the  cube  of 
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the  dimensionality  of  the  matrix.  On  present  day  com¬ 
puters  several  seconds  of  computation  time  may  be 
necessary  to  determine  the  inverse  of  a  l+xl+  matrix# 

This  being  the  case,  one  can  easily  see  that  the  sampling 
rate  may  be  adversely  affected#  In  therefore  behoves 
one  to  make  a  close  study  of  this  matrix  to  see  if  com¬ 
putation  time  can  be  reduced. 

To  illustrate,  let  us  assume  we  are  going  to  track  an 

object  in  space,  receiving  position  information  only.  The 

observable  states  are  range,  bearing,  and  elevation  (R,©, 

§),  and  we  wish  to  determine  best  estimates  of  these 

states  along  with  their  derivatives  (R,  ©,  (J>).  Our  state 

vector  would  then  be  set  up  as  follows: 

R  xl 

R  x2 

©  x3 

x(t)  -  ©  -  x4 

<$>  x5 

|  x6 

The  H  matrix  would  become 

100000 
H(t)  -  001000 

OOQQtO 
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and  hence  the  matrix  (  HPH  +  R  )  would  be  3*3  in¬ 
dimension. 

Recalling  that  P(t)  is  symmetrical,  computational 
reduction  might  be  possible  if  we  can  assume 


E[{X i  ( t )  -£i -° 

-for  L-  i,4 

and  further  that  R(t)  has  non-zero  elements  along  the 
main  diagonal  only.  Making  the  above  assumptions  the  P 
and  R  matrices  would  become 


P(t) 


and  R(t) 


Pn 

p21 

0 

0 

0 

0 

rll 

0 

0 


pl2 

p22 

0 

0 

0 

0 

0 

r22 

0 


0 

0 

P33 

p43 

0 

0 

0 

0 

r33 


0 

0 

p34 

p44 

0 

0 


0 

0 

0 

0 

p55 

p65 


o 

o 

0 

0 

p56 

p66 


We  would  then.'^ind  that 


(  HPH  +  R  ) 


-1 


pll  +  rH 
0 


p33  +  r22 
0 


p55  +  1*33 
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since  the  inverse  of  a  matrix,  having  non-zero  elements 
along  the  main  diagonal  only,  is  found  simply  by  inverting 
the  diagonal  elements. 

Suppose  now  chat  rate  information  is  also  available 
so  that  measurements  of  all  six  states  are  made.  If, 
in  addition  to  the  above  assumptions,  we  can  assume  that 
the  cross— variance  elements  in  the  P  matrix,  involving 
the  even  subscripted  states,  are  also  zero,  then 


where  aii  -  pii  +  rii,  i  =  1,6 

thus  reducing  computational  time  by  at  least  6^/3x2-^  or 


9  times. 

A  further  reduction  might  be  realized  in  the  given 
example  by  the  use  of  three  filters  doing  2x2  matrix 
manipulations  as  opposed  to  one  filter  computing  at  the 
6x6  level. 

Another  consideration  is  the  time  lag  between  input 
and  output.  In  real  time  situations  this  could  be  of  the 
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utmost  importance.  A.  glance  at  the  filter  equations,  I-V, 
indicates  that  if  all  five  equations  are  computed  after  the 
input,  _z(t),  has  been  received,  then  a  considerable  time 
lag  could  ensue  before  the  output  is  generated.  On  the 
other  hand,  if  the  transition  matrix,  ^  ( t,  t— 1 )  ,  is  known 
at  time  t-1,  then  equations  I,  II,  III,  and  IV  can  be 
computed  prior  to  receiving  the  input.  The  time  delay  in 
this  case  would  only  involve  the  time  taken  to  compute 
V,  which  might  be  in  the  order  of  micro-seconds. 

3-2  FILTER  FOR  A  STATIONARY  PROCESS. 

It  is  to  be  noted  that  equations  I  and  II  do  not 
involve  the  observations,  z_(t).  If  the  process,  which  we 
are  trying  to  observe,  is  stationary,  I.ej<^,  H,  Q,  and 
R  are  constant  matrices,  then  It  will  be  found  that  the 
optimal  gain  matrix,  G,  will  stabilize  to  a  constant  value. 
This  matrix  could  then  be  precomputed  (prior  to  any 
observations)  ,  and  the  filter  would  be  reduced  to  the 
relatively  simple  calculations  indicated  by  III,  IV,  and  V. 

A  digital  program,  which  computes  the  optimum 
gain  matrix  for  a  stationary  dynamic  process,  has  been 
written  in  FORTRAN,  and  is  found  in  Appendix  II,  It 


is  completely  general  and  can  handle  any  system  with  up  to 
twelve  states.  It  is  written  as  a  subroutine  to  eliminate 
possible  conflict  in  the  naming  of  variables.  Essentially 
the  subroutine  carries  out  the  iterations  indicated  by  I 
and  II  until  such  time  tnat  no  element  of  the  gain  matrix 
changes  from  its  previous  value  by  an  amount  more  than 
.00001.  If  higher  accuracy  is  desired,  this  constant  can 
be  easily  changed  to  the  degree  of  accuracy  required. 
Further  description  on  the  usage  of  this  program  is  found 
in  the  Appendix. 

3-3  THE  GENERAL  FILTER. 

In  the  general  case,  nonr-stationarity  Is  assumed. 
Appendix  III  contains  two  programs,  the  first  of  which 
performs  the  computations  indicated  by  the  five  filter 
equations  after  each  observation.  The  second  program'^ 4- 
allows  for  the  case  when  the  transition  matrix  is  known 
before  an  observation  is  made  and  hence  reduces  the 
time  lag  (previously  discussed)  between  out  put  and  input. 
Further  discretion  of  the  usage  of  these  programs  Is 
contained  in  the  Appendix. 
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CHAPTER  IV 


ILLUSTRATIVE  EXAMPLES 

t 

This  section  of  the  paper  will  deal  With  two  illustra¬ 
tive  examples.  The  aim  of  the  first  example  is  primarily 
tutorial.  A  thorough  discussion  of  the  model  will  be  given, 
followed  by  the  design  of  an.  optimal  filter.  The  second 
example  will  deal  with  the  track  smoothing  of  an  anti-ship' 
missile.  A  mathematical  model  of  the  missile  in  flight  will 
be  set  up  and  a  filter  designed  for  optimum  track  smooth¬ 
ing. 

EXAMPLE  I: 

/ 

Given  the  transfer  function  for  a  low  pass  filter  in 
Fig.  1,  a)  determine  all  mathematical  model  parameters  , 
and,  b)  design  an  optimal  filter  which  will  give  a  best 
estimate  of  the  states  in  the  filter.  Assume  that  the 

i 

excitation  at  the  input  (u{t)},  and  the  additive  noise 
(v(t))  at  the  output  are  gaussian  and  stationary. 


Fig.  4-1  Cow  pass  filter  of  Example  I. 
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Our  first  step  is  to  convert  Fig..  4-1  to  &  more 


Fig.  4-2  Signal  flow  graph  for  Example  I 
From  Chapter  I,-  equations  (1.8)  through  (1.12),  we  see 
that  x(t)  =  Fx(t)  +  Du(t)  (4.1) 

or 


We  wish  now  to  find  the  solution  to  (4.1)  in  the 
form  given  by  (1.15)  .  To  do  so  we  must  assume  that 
u(t)  is  piece- wise  constant. 

The  and  matrics  will  be  obtained  by  using  signal 
flow  techniques  and  applying  unit  impulses  at  selected 
locations  in  the  diagram  as  illustrated  in  Fig.  4-2.  . 

From  Fig.  4-2  we  see  that 


'■($)  x,(s)" 

61  Sz 
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Taking  the  inverse  of  the  above  and  letting  the  sample 
interval  be  T  we  find 


-T  -ar 

e  -e 


-IT  -T 

le  -ae 


e“T  -e“2T 


zeiT  -  e~r 


(4-2) 
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The  solution  is  now  in  the  form  of  (1.15) 

x(t+r)  =  x(t)  f  A(tfTjt)uCt)  (1.15) 


and  the  mathematical  model  is  shown  in  Fig.  4-3» 


Figure  4-3  Mathematical  model  for  Example  I 

where  H(t)  is  obviously  [l  oj,  since  only  one  state,  xj  , 

observed. 

b)  Design  of  optimal  filter 

We  recall  that  optimum  filtering  is  based  upon  a 

knowledge  of  the  statistic^  of  u(t)  and  v(t),  and  there- 

T"  -r 

foi'e  we  assume  that  E^u(t).u  (t)]  and  E£v(t)'v  (t)Ja: 
part  of  the  problem  statement.  We  now  calculate  the 
covariance  matrix 

Q  (t  +t)  =  A  (t  *  i/(t)]Ar(**r>t) 


(4-4) 


(4.5) 


and  R  (+)  =  £  [v;(t)*  =  n„ 

The  only  remaining  task  is  to  select  an  initial 
covariance  matrix,  P(tG),  for  our  best  estimates  of 
the  initial  states  of  the  model.  The  selection  ■will  depend 
to  a  great  extent  on  a  good  knowledge  of  the  problem. 

In  this  instance  the  best  selection  would  probably  be  the 
main  diagonal  of  Q(t)  (previously  calculated)!  with  the 
off-diagonal  terms  set  to  zero.  The  off-diagonal  terms 
are  zero  initially  since  knowledge  of  xq(tQ)  (at  the  first 
measurement)  will  in  no  way  provide  any  information 
about  x2(t0)  ieXi(tQ)  andXg(t0)  are  uncorrelated  inso¬ 
far  as  the  observer  is  concerned. 

Since  the  system  under  study  is  stationary,  the 
optimum  gain  matrix  may  be  pre-determined.  This 
entails  the  use  of  ■  SUBROUTINE  CONFIL  in  Appendix 
II.  A  sample  period  of  0.1  secs,  is  used.  Using  (4*3) 
and  (4*4)*  we  compute  the  covariance  of  the  states* 
Q(t).  The  element  qjj  in  Q(t)  is  a  measure  of  the 
expected  signal  power.  By  making  r^  (measurement 
noise  power)  equal  to  various  multiples  of  (Including 
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zero)  we  are  able  to  study  the  behaviour  of  the  gain 
matrix  as  a  function  of  noise-.to-signal  power. 

The  curues  in  Figure  4-4  indicate  how  the  optimum 
gain  matrix  elements  behave  as  noise-tc-signal  power  is 
increased.  ‘With  N/S  =  0  we  see  that  <31  is  equal  to 
unity.  We  expect  this  since  there  is  no  measurement 
noise  and  hence  the  best  estimate  of  the  observable 
state  is  the  measurement  z'(t)  itself.  However  as  the 
noise  power  is  increased  the  gain  element  falls  off  and 
the  filter  starts  to  rely  more  on  the  predicted  value 
and  less  on  the  observed  value.  The  matrix  element  G2 
also  falls  off  in  a  similar  fashion,  with  x  2(t)  becoming 
less  dependent  on  z(t)  as  the  relative  noise  power 
increases. 

EXAMPLE  II 

Problem  Statement  -  It  is  known  that  the  enemy's 

/ 

main  anti~ship  weapon  is  an  air-launche.d  missile  which  is 
normally  launched  at  a  distance  of  250  to  300  miles  from 
target.  After  launch,  it  climbs  to  an  altitude  of  40,000 
ft.,  attains  a  speed  of  approximately .  1000  mph,  and 
maintains  this  speed  by  use  of  a  sustainer  motor.  When 
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•within  25  miles  of  the  expected  location  of  target  a 
search  device  is  switched  on  which  pinpoints  the  target 
and  enables  the  missile  to  guide  itself  to  target.  A 
typical  friendly  ship  at  sea  is  fitted  with  a  mono-pulse 
search  radar  system  with  a  scan  rate  of  10  scans/ min. 
The  ship  is  fitted  with  a  digital  computer  and  target 
information  is  available  to  the  computer.  It  is  known 
..that  the  probability  distribution  function  of  the 
accumulated  error  (in  both  x  and  y)  is  approximately, 
normal  with  zero  mean  and  2  mile  standard  deviation. 

.  Data  accumulated  on  similar  missiles  indicate  that, 
due  to  erratic  thrust  developed  by  the  missile  motor  and 
average  atmospheric  turbulence,  the  velocity  of  the 
missile  varies  in  a  random  fashion  (approximately 
gaussian) .  The  standard  deviation' of  this  randomness 
is  about  Z%  of  mean  velocity.  __ 

Design. a  filter  which  will  determine  best  estimates 
of  the  missile’s  position  .and  velocity.  Assume  that 
attack  is  equi-likely  from  all  directions. 

SOLUTION 

Our  first  step  is  to  set  up  a  mathematical  model 
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which  describes  the  dynamics  of  the  missile,  in  the  form 


x  (±)  =  F (t)  X  (t)  t  DC-t)u(t) 


Considering  only  the  component  in  the  X  direction  we 


obtain 

X,  (t) 

*“  r-* 


> 


where  Xi  is  the  position  on  the  x  axis 

x^  is  the  velocity  component  in  the  x  direction 
A  similar  vector  equation  would  describe  the 
dynamics  in  the  y  direction. 

As  we  have  seen  earlier  the  solution  to  the  above 
equation  Is 

X(t)  t,)  +  A  ux(t)  (4-7) 


assuming  ux(t)  Is  piecewise  constant. 

We  now  must  determine  $  ( t^ ,  t  j, ) . 

.From  (4*4)  we  obtain  the  model  shown  in  Fig  3*5* »  & 
simple  1/s1  plant. 


34 


JL 

_L 

X.  . 

"  f 

5 

S 

Fig  '3.5  Model  for  missle  in  Example  II. 
The  parameters  ^  and  A  are 

$(**;*»)  = 
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o  / 


Since  we  intend  to  sample  at  a  rate  of  10  times/ 

min.  (scan  rate  of  radar)  then  tg  -,.tj  «  6  seconds. 

or  T  »  _ 1^  hrS. 
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Our  next  step  is  to  „ find  Q(t)  \ 


X  ^ 

Q  (t-)  =  £(t,t  -r)  E  [ujlt-y)- uJtt-Tj]  Ar 
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The  problem  statement  specified  however  that  ^ 
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Summarising  our  model  we  have  (Fig.  4.6) 


Fig  4*6  Model  for  Missis  of  Example  II 

We  are  now  ready  to  set  up  the  filter.  One 
requirement  is  the  initial  covariance  matrix  P(o).  Since 
the  missile  can  approach  us  from  any  direction  with 
equal  probability  our. best  estimate  for  the  inital  states 
is  zero  for  both  position  and  velocity,  ie 
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We  therefore  find  that 


P(o)  -  E[2(»)-h0 j 


s, 

a  \ 


To  determine  this  we  must  have  a  knowledge  about 
the  detection  capability  of  the  radar.  A  standard  search 
radar  looking  for  an  object  at  40,000  ft.  elevation  may 
have  an  initial  detection  density  function,  as  shown. in 
Figure  4*7 ,  in  all  directions  in  the  x-y  plane. 


0  <r=2  00 


Range  (miles) 


FiG.  4.7  Probability  density  function  for  inital 
detection 

Since  the  missile  may  approach  from  any  direction 
let  us  assume  that  the  standard  deviation  when  con¬ 
s' iering  the  x  direction  only-  is.  100  miles.  Similarly  let 
us  assume  the  velocity  component  on  the  x  axis  has  a 
standard  deviation  of  500  mph,  ,  We  then  find  that 


P(O) 
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Our  last  remaining  task  is  to  select  a  value  for  Vx^(G) 
to  provide  numbers  for  the  Q(o)  matrix.  Our  best 
initial  estimate  of  this  quantity  would  be  the  pgg  element 


in  F(o)  . 
Hence 


V  ■ 


(o)  =  ( SooY 


Jlu  0 


Subsequent  values  of  V„(t)  could  be  calculated  by 


A 


using  x2(t) 


We  are  now  prepared  to  write  a  program  for  the 


optimal  filter. 


Since  the  transition  matrix  is  a  constant  but  .  Q 
is  variable  (being  a  function  of  Vx)  the  2nd  program  of 
Appendix  III  would  apply. 

We  now  summarize  taking  into  account  the  y  com¬ 
ponent  of  direction.  The  flow  chart  for  the  filter 


calculations  is  shown  in  Appendix  III. 
n  (number  of  states)  “  1+ 

P  (number  of  observables)  »  2 
T  (sample  interval)  =  1/600  hrs. 


39 


x(  t) 

x1(  t) 

x(  t) 

*2(t> 

y(t) 

x3(  t) 

y(t) 

x4(t) 

£(t) 


x1(t) 

x3(t) 


*l(t)  +  /x(t) 

x3(t)  ♦:  y(t) 


r 


/ 


H(t) 


10  0  0 
0  0  10 


$(t+T,t)  = 


1  T  0  0 

0  10  0 
0  0  1  T 

0  0  0  1 


°*22(°)  “  <144(0) 


“  ( .02  x  50 0)2 


100 


40 


BIBLIOGRAPHY 


1.  Kalman,  R.  E.;  "A  New  Approach  to  Linear  Filtering 
and  Prediction  Problems";  .Journal  of  Basic  Engineer¬ 
ing,  March,  i960. 

2.  Kalman,  R.  E0;  "New  Methods  and  Results  in  Linear 
Prediction  and  Filtering  Theory";  Rias  Technical 
Report  6l-l,  1961. 

3.  Magill,  T.  D0;  "Optimal  Adaptive  Estimation  of 
Sampled  Stochastie  Processes";  Stanford  Electronics 
Laboratories  Technical  Report  No.  6302-2,  December 
1963. 

4.  Schmidt,  S.  F.;  "The  Application  of  State  Space 
Methods  to  Navigation  Problems";  Philco  Technical 
Report  No.  4»  ^July  1964. 


43 


APPENDIX  I 


LINEAR  DIFFERENTIAL  EQUATIONS 
Let  us  consider  the  solution  of  a  1st  order  dif¬ 
ferential  equation  by  the  use  of  an  integrating  factor 
(p),  to  make  an  exact  differential. 

Given: 


dx 

dt 


+  ax 


u 


Find: 


x(t)  =  f(xQ,  U,  a,  t) 


Solution: 

Multiply  (1)  by  p(t)  and  attempt  to  make  the 
resulting  equation  an  exact  differential.  We  have 

dx 

P-$r  +  pax  -  pu 


or 


(2) 


dp 

x  +  aPx  "  Pu 


and 

ft  (px)  =  (dl  - ap)  x  +  pu  (3) 

Considering  the  homogeneous  part  of  the  problem 
(i.e.,  u  =  0),  we  can  make  (3)  an  exact  differential 
by  setting 


(4) 


4:5 


represents  the  particular  solution  or  convolution  integral* 
Now  let  us  consider  a  set  of  n  of  these  1st  order 
differential  equations 
do-:. 

dF  *  f^*)'  u)  i  -  1,  ....  .7 

or 

x  =  Foe  +  Du  (11) 


where  x  and  u  are  1  x  n  column  vectors  and  P  and  D 
are  n  x  n  matrices.  Multiplying  (11)  by  the  integrating 
factor  'p(t),  (a  row  vector)  we  obtain,  after  some 
manipulation 


d 

dt 


(  V  x  )  .  * 


(12) 


As  before  we  assume  a  solution  for  the  adjoint  variable, 
p(t)  as 

_  _  -Ft 

P  -  P0e 

.Substituting  into  (12)  and  integrating  gives 
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o  L 


-Ft 


x  -  x  * 

—  — o 


J  e"F^Du(r)  dr 


(13) 


Ft 

Multiplying  both  sides  by  e  gives 


Ft 

x  =  e  x  ,  \  e 

—  — O  T 

tfo 


f  f  eF(t  ~r)  Du  (r)  dr 


(14) 


Often  the  above  is  expressed  in  the  form 

xttj).  =  $(tj,  tQ)  x(tQ)+  A(tlf  t0)  u(t0)  (15) 

V 

Where  (£ is  the  transition  matrix  or  fundamental  matrix, 
and  A  represents  the  distribution  matrix  with  u(*T) 
held  constant  at  its  value  u(tQ). 
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APPENDIX  II 


THE  STATIONARY  FILTER 

PURPOSE: 

The  attached  program  can  be  used  to  determine  the 
stabilized  optimum  gain  matrix  for  the  estimation  of  the 
states  of  a  stationary  process, 

USAGE : 

1.  .Calling  Sequence 

CALL  CONFIL  (P,  Q,  R,  TR,  H,  KN,  KP,  KER, 

G) 


2.  Arguments 

a.  P  -  .initial  covariance  matrix  of  system  states  - 

dimension  (12  x  12) 

b.  Q  -  covariance  matrix  .of  states  due  to  gaussian 

excitation  -  dimension  (12  x  12) 

c.  R  -  covariance  matrix  of  measurement  noise  - 

dimension  (12  x  12) 

d.  TR  -  transition  matrix  of  pro, cess  - 

dimension  (12  x  12) 

e.  H  -  matrix  which  defines  the  observable  states  - 

dimension  (12  x  12) 

f.  KN  -  number  of  states  in  the  system - 

dimension  (scaler) 
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g.  KP  -  number  of  observable  states  - 

dimension  (scaler) 

h.  KER  -  error  indication  (=  2  implies  the  inverse 

of  a  matrix  could  not  be  obtained)  — 
dimension  (1) 

i.  G  —  optimal  gain  matrix  - 

dimension  (12  x  12) 

3.  Accuracy:  .see  Chapter  III 

2+.  Equipment  Configuration:  CDC  l601f  with  FORTRAN 

t 

60. 

5*  Cautions  to  Users 

a.  All  arguments  in  the  main  program  must  be 
dimensioned  the  same  as  those  in  CONFIL. 

b.  The  main  program  should  contain  an  ERROR  TEST 
( see  2h. ) 

6.  Flow  chart  showing  Typical  Usage  (see  Fig.  11-.1)  . 


Fig.  11-1  Flow  diagram  showing  typical  usage 
of  CONFIL 
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APPENDIX  II  (CONTINUED) 


APPENDIX  II  II-  SUBROUTINE  CONFIL  (CONTINUED) 


END  STABILIZATION  TEST 


APPENDIX  II  II-  SUBROUTINE  CONFIL  (CONTINUED) 
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APPENDIX  II  II-  •  SUBROUTINE  CONFIL  (CONTINUED) 


APPENDIX  III 


GENERAL  FILTERS 

PREAMBLE:  Two  programmes,  in  subroutine  form,  are 

contained  in  this  appendix.  ..The  first  programme  carries 
out  all  filter  computations  after  the  input  has  been 
received  at  each  sample  instant.  The  time  lag  between 
output  and  input  will  therefore  depend  on  the  time  taken 
for  these  computations.  The  second,  programme  is 
designed  for  use  when  certain  parameters  are  known  (or 
at  least  a  very  good  estimate  of  these  parameters  can 
be  made)  prior  to  receiving  the  input  at  each  sample 
Instant.  A  prior  knowledge  of  these  parameters,  namely 
TR(  tpt—l) ,  H(t),  R(t),  Q(t)  (defined  below),  enables  a 
considerable  reduction  of  the  above-mentioned  time  lag. 
between  output  and  input.  This  is  achieved  by  perform¬ 
ing  most  of  the  filter  computations  prior  to  receiving 
the  input. 

A.  SUBROUTINE  BESTX 

PURPOSE:  This  subroutine  will  provide  an  optimum 

estimate  of  the  state  vector  for  any  sampled-data 
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linear  dynamic  process  (with  twelve  or  less  states)  if 
booh  the  process  random  excitation  and  corruptive 
measurement  noise  vectors  have  gaussian  distributions. 
USAGE: 

1.  Calling  Sequence  * 

CALL  BESTX  (  F,  Qa  R,  TR,  H,  KN,  KP,  KER, 
Gs  XP,  Es  X) 

2.  Arguments 


a.  P  -  initial  covariance  matrix  of  system  states  - 

dimension  (12  x  12) 

b.  Q  -  covariance  matrix  of  states  due  to  gaussian 

excitation  - 
dimension  (12  x  12) 

c.  R  -  covariance  matrix  of  measurement  noise  - 

dimension  (12  x  12) 

d*.  TR  -  transition  matrix  of  process  - 
dimension  (12  x  12) 

e«  H  -  matrix  which  defines  the  observable  states  - 
dimension  (12  x  12) 

f «  KN  -  number  of  states  in  the  system  - 
dimension  (scaler) 

g.  KP  -  number  of  observables  states  - 

dimension  (scaler) 

h.  KER  -  error  indication  (-2  implies  the  inverse  of  a 

matrix  could  not  be  obtained)  - 
dimension  (1) 

i.  G  -  optimal  gain  matrix  - 

dimension  (12  x  12) 


56 


XP 


j  * 


-  predicted  estimate  of  process  state  vector  — 
dimension 


k.  3  -  observation  vector  - 
dimension  (12  x  12) 

1.  X  -  optimal  estimate  of  process  state  vector 
(generated  by  filter)  - 
dimension  (12  x  12) 


3*  Equipment  Configuration;  CDC  1604  with  FORTRAN 


60. 


4*  Cautions  to  User; 

a.  All  arguments  in  the  main  program  must  be 
dimensioned  the  same  as  those  in  BESTX 

b.  The  main  program  should  contain  an  ERROR  TEST 
( see-  2h, ) 

c.  See  attached  flow  chart  depicting  typical  usage. 
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APPENDIX  III  III-  SUBROUTINE  BESTX ( CONTI NUED ) 


(JN3 

00  VO 0  -  0  5 

OQVOO  NJjni3y 

OEVOO 


SUBROUTINE  BESTXtCONT INUED) 


B.  SUBROUTINE  GANDPR 

PURPOSE:  This  subroutine  will  provide,  for  each 

sample  instant,  an  optimum  gain  matrix,  and,  the  pre¬ 
dicted  values  of  the  process  and  the  observation  state 
vectors,  if  the  process  is  linear  and  the  random 
excitation  and  corruptive  measurement  noise  vectors 
have  gaussian  distributions. 

USAGE: 


1.  Calling  Sequence 

CALL  GANDPR  (P,  Q,  R,  TR,  H,(X,  XP,  2P, 
KN,  KF,  KER,  G) 


2.  Arguments 

a.  P  -  initial  covariance  matrix  of  system  states  - 

dimension  (12  x  12) 

b.  Q-  covariance  matrix  of  states  due  to  gaussian 

excitation  -  dimension  (12  x  12) 

c.  R  -  covariance  matrix  of  measurement  noise  - 

dimension  (12  x  12) 

d.  TR  -  transition  matrix  of  process  - 

dimension  (12  x  12) 


e.  H  -  matrix  which  defines  the  observable  states  - 
dimension  (12  x  12) 
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f .  X  -  optimal  estimate  of  process  state  vector  — 

dimension  (12  x  12) 

g.  XP  -  predicted  estimate  of  process  state  vector  - 

dimension  (12  x  12) 

h.  SP  -  predicted  estimate  of  observable  state 

vector  - 

dimension  (12  x  12) 

i.  KN  -  number  of  states  in  the  system  — 

dimension  (scaler) 

j.  KP  -  number  of  observable  states  - 

dimension  (scaler) 

lc.  KER  -  error  indication  (=>  2  implies  the  inverse 
of  a  matrix  could  not  be  obtained)  - 
dimension  (1) 

1.  G  -  optimal  gain  matrix  - 
dimension  (12  x  12) 

3.  Equipment  Configurations  CDC  1604  with  FORTRAN 


4.  Cautions  to  User: 

a.  All  arguments  in  the  main  program  must  be 
dimensioned  the  same  as  thoge  in  GANDpR. 

b.  The  main  program  should  contain  an  ERROR  TEST 
( see  2h„ ) 

c.  See  attached  flow  chart  depicting  typical  usage. 


111-2  Flow  diagram  showing 

typical  usage  of  GANDPR 
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.APPENDIX  III  III-  SUBROUTINE  GANDPR(CONTINUED) 


SUBROUTINE  PROD  t  A,  B , N  ,M,L  , C  I 
DIMENSION  A<12»12), 6(12*12), G(12, 12) 
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APPENDIX  III  III-  SUBROUTINE  GANDPR ( CONT I NUED ) 
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SUBROUTINE  GANDPR(CONTINUED) 
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